% QDMBzFromBu.m -- Auxiliary MATLAB script for computing the Bz field map
%    from the B111 (Bu) field map in the Fourier domain. This script is called
%    by QDMdataprocessing.m
%
%  References: - R.R. Fu, E.A. Lima, M.W.R. Volk, R. Trubko (2020) "High-sensitivity
%  moment magnetometry with the quantum diamond microscope,"  Geochemistry,
%  Geophysics, Geosystems.
%              - E. A. Lima and B. P. Weiss (2009) "Obtaining vector magnetic field
%  maps from single-component measurements of geological samples," Journal
%  of Geophysical Research - Solid Earth, 114, B06102.
%
% -----------------------------------------------------------------------------------
%                Matlab code by Eduardo A. Lima and Roger R. Fu
%      Copyright (C) 2017-2020 MIT Paleomagnetism Lab, Harvard Paleomagnetics Lab
% -----------------------------------------------------------------------------------

function [Bz] = QDMBzFromBu(Bu,fs,u)
%   ---------------------------------------------------------------------------------
%   Bu     -> u component of the magnetic field measured in a regular planar grid
%   fs     -> Sampling frequency in 1/m
%   u      -> unit vector representing the direction of the measured field component
%   ---------------------------------------------------------------------------------

SHOWGRAPHS=0;      % Set this value to 1 to show frequency response plots, or to 0 otherwise.

[SIZEx, SIZEy] = size(Bu);
N1=SIZEx;
N2=SIZEy;

f1 = [0:N1/2-1 -(N1/2):-1]*fs/N1;         % These freq. coordinates match the fft algorithm
f2 = [0:N2/2-1 -(N2/2):-1]*fs/N2;

[F2,F1] = meshgrid(f2+1e-30,f1+1e-30);    % Add an epsilon to avoid singularity at the origin

ky = 2*pi*F1;
kx = 2*pi*F2;

k=sqrt(kx.^2+ky.^2);


etz=k./(u(3)*k - u(2)*1i*ky - u(1)*1i*kx);    % Calculate the filter frequency response
                                              % associated with the z component

e = fft2(Bu, N1, N2);

Bz=ifft2(e.*etz,N1,N2,'symmetric');           % Compute Bz map from Bu map



